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Abstract 

We investigate the effect of random columnar disorder on the superconducting phase transition 
of a type-II superconductor in zero applied magnetic field using numerical simulations of three 
dimensional XY and vortex loop models. We consider both an unscreened model, in which the 
bare magnetic penetration length is approximated as infinite, and a strongly screened model, 
in which the magnetic penetration length is of order the vortex core radius. We consider both 
equilibrium and dynamic critical exponents. We show that, as in the disorder free case, the 
equilibrium transitions of the unscreened and strongly screened models lie in the same universality 
class, however scaling is now anisotropic. We find for the correlation length exponent v = 1.2ib0.1, 
and for the anisotropy exponent C, = 1.3 ib 0.1. We find different dynamic critical exponents for 
the unscreened and strongly screened models. 

PACS numbers: 
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I. INTRODUCTION 



The discovery of high temperature superconductors, in which thermal fluctuations are 
important and mean field theory can no longer be applied, has led to a resurgence of inter- 
est in phase transitions and critical phenomena in type-II superconductors.— Of particular 
interest has been understanding the effects of random quenched disorder on the nature of 
the ordered phases and the universality of the phase transitions. For the high temperature 
superconductors, this quenched disorder can take many forms: random point defects due 
to oxygen vacancies, planar twin grain boundaries, and columnar defects introduced by ion 
irradiation. 

Most of the work in this area has focused on the case of a finite applied magnetic field, 
where one seeks to understand how the randomness distorts or destroys the Abrikosov lattice 
of magnetic field induced vortex lines that forms in a pure system. Columnar defects^i^ 
have received considerable attention, as they are particularly effective in pinning vortex lines 
and reducing flux flow resistance. In contrast, in this paper we will focus on the effect of 
columnar defects on the superconducting phase transition in zero applied magnetic field. We 
expect this case to be interesting for the following reason. A generalized Harris criterion^'^ 
argues that disorder will be a relevant perturbation, and change the universality class of a 
phase transition, whenever 2 — d*v > 0, where d* is the number of dimensions in which the 
system is disordered, and u is the usual correlation length critical exponent. For a disorder 
free superconductor, the transition in zero applied magnetic field is in the universality class 
of the three dimensional (3D) XY model,— for which u > 2/3. For random point disorder, 
d* = 3, so 2 — (i*z/ < 0, and the generalized Harris criterion argues that the universality of the 
transition remains unchanged. For columnar disorder, however, d* = 2, and so 2 — d*^ > 0. 
Columnar defects should therefore cross the zero field transition over to a new universality 
class. Stability^ of this new disordered fixed point with respect to the generalized Harris 
criterion implies that it should have a new value u > 1. In our equilibrium simulations we 
indeed find behavior consistent with this, and we obtain a value for the correlation length 
exponent u = 1.2 ± 0.1. Moreover, we find scaling is now anistropic and we find the value of 
the anisotropy exponent to be C = 1.3 ±0.1. Experimental measurement of these exponents 
would therefore provide a precision test of the theoretical model. 

The model we study also has application to the T = superconductor to insulator 
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quantum phase transition in two dimensional thin films with random substrate disorder , 
and to the Mott transition for bosons in 2D optical lattices with the addition of random 
scattered laser intensity.— In these cases, the two dimensional quantum problem can be 
mapped onto a corresponding three dimensional classical problem with the same symmetries 
as the one we study here.— 

To study the effect of columnar disorder on the zero field transition of a type-II super- 
conductor, we will consider two different limits. The first is the limit of an "unscreened" 
superconductori^^ in which magnetic field fluctuations are frozen out, corresponding to the 
approximation of an infinite bare magnetic penetration length, Aq ^ c>o. Here, vortex line 
segments have long ranged Coulombic-like interactions. For the extreme type-II high tem- 
perature superconductors, for which k = Aq/^o ~ 100, where ,^0 is the bare coherence length 
that also sets the radius of a vortex line core, the unscreened model should give a good 
description except in an extremely narrow temperature window about the transition.— 

The second limit is that of a strongly screened superconductor,i^*i^ corresponding to the 
case Ao ~ ^o- In this case, vortex line segments have short range interactions. This de- 
scription should also become valid extremely close to the transition when the diverging 
correlation length ^ becomes comparable to the renormalized magnetic penetration length 
•^5 ^ and magnetic field fluctuations on such large length scales must be included in 
determining the true critical behavior. This region near Tc may, however, be too small to 
observe in practice.— 

As in the disorder free duality transformationi2*iSiil establishes that these two 

limits lie in the same universality class as regards equilibrium critical behavior. They may be 
different, however, for dynamic critical behavior.— In this work we carry out detailed Monte 
Carlo (MC) simulations of the XY model for the unscreened superconductor to determine 
the equilibrium critical exponents, and we demonstrate the presence of anisotropic scaling; 
by duality, these exponents also apply to the strongly screened case. Then, using simple local 
Monte Carlo dynamics, we compute the dynamic critical exponent for both the unscreened 
XY model, and for the strongly screened vortex loop model. We find that the dynamic 
exponent is different for these two limits. 

The remainder of this paper is organzied as follows. In Section II we describe the XY 
model and the loop model for the unscreened and strongly screened limits, respectively. The 
duality transformation between the two is given in Appendix A. In Section III we discuss 
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the equilibrium critical behavior of the XY model, presenting our finite size scaling analysis, 
defining the observables we measure, and giving the numerical results of our simulations. 
In Section IV we discuss the dynamic critical behavior of the XY and loop models, within 
a simple local Monte Carlo dynamics. We define the observables we measure and give our 
numerical results. In Section V we give our discussion and conclusions. 

II. MODELS 
A. XY Model 

To model the effects of thermal fluctuations in a type-II superconductor, we start with the 
commonly used 3D XY model.— This models the phase fluctuations of the superconducting 
order parameter in the "unscreened" limit where magnetic fleld fluctuations are frozen out, 
corresponding to the approximation of an inflnite bare magnetic penetration length, Aq — > oo. 
For zero applied magnetic fleld we have. 



Here 6i represents the phase angle of the complex superconducting order parameter on node 
z of a periodic cubic grid of N = L x L x Lz sites, with periodic boundary conditions in all 
directions. The sum is over all nearest neighbor bonds of the grid, with fi = x,y,z, 

and the cosine term represents the kinetic energy of fluctuating supercurrents. The short 
length cutoff of the discrete grid models the bare vortex core size ^o- 

In a pure system, the couplings Jj^ are all equal, except for a possible variation with 
bond direction /i. Here, we take the Jj^ randomly distributed in order to model quenched 
random columnar defects. For the work reported on here, with columnar defects aligned 
parallel to the z axis, we have chosen the following distribution: in the z direction, we take 
all Jiz = 1; in the xy plane, we take J,^, fi = x,y, distributed equally likely with the two 
values 0.1 and 1.9, keeping the Jj^ translationally invariant along the z axis so as to model 
columnar disorder. Note that the random Jj^ introduce no frustration into the system; in 
the ground state all the 6i are equal. The variations in the Jj^ result in spatially random 
pinning energies for vortex loop excitations of the phase angles 6i. We have chosen the above 
bimodal distribution for Jj^ to give strong pinning energies (for flxed average Jj^), so as to 
be able to approach the asymptotic scaling limit with reasonable size systems. 




(1) 
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Although we will simulate the Hamiltonian of Eq. (^J) using periodic boundary conditions 
on the phase angles 6'j, it is useful to consider a more general fixed twist houndaij condition, 

= + , (2) 

where is a fixed (non-fluctuating) total twist in the phase angle applied across the system 
in direction fi. Periodic boundary conditions correspond to the twist A^ = 0. Transforming 
to new variables, 

e: = 9,-iA,/L^)r,-fi , (3) 
the Hamiltonian of Eq. becomes, 

HxYie',; \] = -Yl - - ^JL,) , (4) 

where the 9[ obey periodic boundary conditions, O'i^L^n = 0[. Using the fact that the cosine 
is periodic in 27r, the partition function integrals over 6[ can be taken over the interval 
9[ G [0, 27r), as were the integrals over the original phase angles 9i. Considering how the free 
energy varies with the twist A^ will be useful later for discussing phase coherence in the 
model. 



B. Loop Model 

Although we carry out our equilibrium simulations directly in terms of the XY model 
of Eq.(^, we also consider a different formulation of the model. If instead of the cosine 
interaction of Eq.(^, one uses the periodic Gaussian interaction of Villain,— then a standard 
duality transformation^iiSiil (see Appendix A) maps the XY model, TYxy/^) onto a model 
of sterically interacting loops, Tiioop/T, where. 

The rij^ are integer valued variables on the bonds (i, yu) and satisfy a divergenceless con- 
straint, 

^ Km - = . (6) 

The thus form connected paths through the system that must eventually close upon 
themselves. The couplings gi^ of Eq.((Hl) are related to the couplings of the XY model by, 

Qii^/T = T/Ji^, , (7) 



where the temperature scale of the loop model, T, is inverted with respect to the temperature 
scale of the XY model, T. 

While the loop model of Eq.(jn)) is exactly dual to the XY model of an unscreened super- 
conductor, taking it on its own with T as the physical temperature, we can give T^ioop the 



following different physical interpretation.^*i^ We can regard the divergenceless variables n. 



as the vortex loops of a strongly screened superconductor with Aq ~ ^o- The short ranged 
vortex line interaction of this case is then modeled by the simple onsite repulsion of Tiioop- 
Further details of this analogy may be found in RefjlSl. If we regard each site of our numer- 
ical grid as representing a columnar pin, the random Qi^ in the xy plane can be thought of 
as modeling the random distances between such pins, and hence giving the random ener- 
gies associated with a vortex loop segment hopping from one pin to another. This duality 
between Tiioop and TCxy thus implies that the unscreened and the screened superconductor 
models fall in the same equilibrium universality class, just as is the case for the disorder free 
model.-^ 



III. EQUILIBRIUM CRITICAL BEHAVIOR 

In this section we report on our equilibrium XY model simulations. To extract critical 
exponents, we use the method of finite size scaling. We first, therefore, discuss this method. 

A. Finite Size Scaling 

Because the columnar disorder singles out the special direction z, we must allow for the 
possibility that scaling will be anisotropic. If ^ denotes the correlation length in the xy 
plane, then anisotropic scaling assumes that, as T — > Tc and ^ diverges, the correlation 
length along the z axis diverges as, 

~ (8) 

where ( is the anisotropy exponent. 

Consider now an observable O whose scaling dimension is zero. As a function of reduced 
temperature t = {T — T^ /T^. and system size Lx L x L^, we expect the scaling relationship, 

0{T, L, L,) = ditb'/", L/b, Ljb<) , (9) 
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where b is an arbitrary length rescahng factor, O is the scahng function, and u is the usual 
correlation length exponent, 

e ~ t-v (10) 

Taking b = L in Eq,® then gives, 

O(T,L,L,) = a(tLl/^l,L,/L0 . (11) 

For the case of isotropic scaling, with C = 1? choosing a constant aspect ratio = jL reduces 
the right hand side of Eq.()lH) to a function of the single scaling variable tL^/'^ . Measuring O 
vs. T for systems with varying L but fixed L^/ L is then sufficient to determine the exponent 
V. However when C 7^ I5 and its value is unknown, it becomes necessary to consider systems 
with varying aspect ratio L^/L, greatly increasing the complexity of the computation. 

To deal with this case we take the following approach, originally used to study the phase 
transition in the quantum Ising spin glass. Assume that the observable 0{T, L, Lz) when 
viewed as a function of L^, for fixed T and L, has a maximum at a particular value L^max- 
Because of the scaling law Eq.(fTT|). this value L^max must occur when 

L,^,JL^ = ^{tL'/'') , (12) 

where 7 is a scaling function of the single variable tL^/'^ . We then define, 

0^ax(T,L) = 0(T,L,L,^aJ = a(^Ll/^l,7(^Ll/'^)) = O^^tL^/'') . (13) 

Plotting Ojnaxl^) L) vs. T for different values of L, the curves will intersect at the common 
point T = Tc (i.e. t = 0). The slopes of these curves at then determine the exponent u. In 
practice, we will determine the values of Tc and the exponent u by the following approach.— 
Close to Tc (i.e for small t) we can expand the scaling function Omax as a polynomial for 
small values of its argument tL^/'^, 

d^^UtL'^n -ao + ai[(T - n)/T,]L'/'' + ([(T - T,)/T,]L^/-f + . . . (14) 

We then fit the data for (9max(T, L) to the above form using Tc, z/, oq, ai, 02 ... as free fitting 
parameters. Varying the system sizes L and temperature window |T — Td of the data used in 
the fit, as well as varying the order of the above polynomial expansion, will give confidence 
on the significance of the fit. 
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Having obtained the value of Tc, plotting I/j;max(^c) vs. L determines the anisotropy 
exponent ( by Eq. fT^ . 

i:.max(r,)~L^ . (15) 

Knowing Tc, u and (, plotting Oraa.xiT,L) vs. tL^^'^ and 0{Tc, L, L^) vs. L^/L'^ should 
collapse the respective data to a single scaling curve. 



B. Observables 



To carry out the scaling analysis outlined in the previous section, we now have to deter- 
mine appropriate observables to measure. 

For the 3D XY model of Eq.([Q), we expect below a non- vanishing order parameter, 
ijj = (l/N) Xli^*^'- We define the real part of tfj as. 



ij^cos^, , (16) 

and construct its Binder ratio— 



^ N 



g{T,L,L,) = l 



^(tLl/^L,/LO . (17) 



_3(M2)2_ 

Because the scaling dimension of M cancels in taking the ratio above, the Binder ratio g 
has scaling dimension zero, and so has the scaling form of Eq. dlip . In the above, (...) 
denotes the usual thermal average, while [. . .] denotes the average over different realizations 
of the columnar disorder. In the denominator of Eq. lfTTj) . the square of the expectation 
value is evaluated using two replicas with identical disorder, indexed by a and b, (M^)^ = 
((M'')2)((M^)2), in order to avoid bias.^^ 

Another observable we have measured is obtained by considering the dependence of the 
total free energy on the total applied twist across the system.— Sensitivity to boundary 
conditions, in this case specified by the twist in Eq.Q, is one of the signatures of an 
ordered phase. The XY model is therefore phase coherent when the total free energy JF 
varies with twist A^. JF is computed from a partition function sum over the 6'- using the 
Hamiltonian T-Cxvl^fi] of Eq.(j3]). A convenient measure of the dependence of JF on A^ is 
obtained by looking at the curvature of JF(A^) at its minimum. In Appendix A we show 
that this minimum always occurs at A^ = 0. We therefore consider. 



dAl I 9A, 
8 



2\ 



(18) 



where T^xy is that of Eq.Q, and the averages on the right hand side are taken in the 
ensemble with = 0. 

Since the total free energy T and the total twist are both scale invariant quantities, 
then S^JF/SA^ has scaling dimension zero. These derivatives are usually defined in terms 
of the helicity modulus^^ T^, which is the derivative of the free energy density with respect 
to the twist per length. We have in three dimensions, 

- ^T. , (19) 



LzT, . (20) 



Averaging over the disorder, we find that, for fixed T and L, {L"^ /Lz)[Tz] decreases mono- 
tonically as increases, while L^[Ta.] increases monotonically as increases. In order to 
have an observable which has a maximum as a function of Lz, we therefore consider the 
product, 

L2[T,TJ =n(^Ll/^L,/LO , (21) 
which has the same scaling form as Eq. ljllj) . 



C. Monte Carlo Methods and Error Estimation 



In order to achieve accurate results, averaging over many disorder realizations for many 
different aspect ratios, Lz/L, it is essential to have an efficient simulation algorithm. For 
equilibrium simulation of the 3D XY model, the lack of frustration allows us to use the 
Wolff^^ cluster algorithm to avoid critical slowing down. We typically use 100 Wolff sweeps 
to approach equilibrium, followed by 200 Wolff sweeps to compute averages; one Wolff sweep 
is defined as building clusters until each phase angle 9i has been updated once on average. 
Between 3000 and 5000 different realizations of the random disorder are averaged over near 
the critical point, with fewer realizations used away from the critical point. A test of the 
equilibration of our simulations is shown in Fig. Q where we see that the above simulation 
lengths are sufficient. 

To estimate the statistical error in our results we use the following method. For our raw 
data, our average is just the average over the individual values obtained in independent 
realizations of the random disorder. Our estimated error is determined from the standard 
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FIG. 1: Binder ratio maximum for the 3D XY model, fl'max at T = 2.05 ~ Tc, computed with N^q 
equihbration sweeps fohowed by A'eq sweeps to compute averages. A'gq = 100 is sufficient for good 
equihbration for ah sizes L. 

deviation a of these A^^^ independent values, error = a /^/Wd- To estimate the statistical error 
in the fitting parameters of our finite size scaling analysis, we take the following approach. 
From our original data set we construct many (typically 1000) fictitious data sets by adding 
to each data point a random Gaussian variable with zero mean, and standard deviation 
equal to the estimated statistical error of the data point. We then fit each of the fictitious 
data sets. The standard deviation of the values of the resulting fitting parameters then gives 
our estimate of the statistical error in the fitting parameter. 

Harder to estimate are the possible systematic errors in our results. Here we rely on 
varying parameters of our analysis, such as the order of a polynomial fit, or the system sizes 
L used in the fit, in order to get a feeling for the likely accuracy of our results. 

D. Results 

We now present our results from simulations of the XY model of Eq.(^. In Fig. |21 we 
plot our results for the Binder ratio of Eq. lfTTj) . g{T, L, L^) vs. L^, for sizes L = 8, 12, 20, at 
the fixed temperature T = 2.05. We see that for each L, g[T, L, Lz) has a clear maximum 
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FIG. 2: Binder ratio, g{T, L, L^) vs. for several values L at the fixed T = 2.05. Solid curves 
are cubic polynomial fits in InL^. 

at a particular Lzma,x- Note that the maximum values of these curves appear to be equal for 
the different values of L. From Eq. (jl3|) we therefore infer that the temperature T = 2.05 
is approximately the critical temperature Tc. To determine the precise values of L^max and 
the maximum values gmax(T,L) = g(T, L, L^max), we fit the data for each L to a cubic 
polynomial in InL^ (these are the solid curves in Fig. 12)). The L^max obtained this way are 
not, in general, integer values. We have also done such fits using a quadratic polynomial 
in InL^; the difference in values obtained from the cubic vs. the quadratic fit provides 
our estimate of the systematic error of this procedure. We find that for QmaxiT, L) this 
systematic error is always smaller than the estimated statistical error of the cubic fits; for 
Lzmax the systematic error is bigger. This reflects the simple fact that g{T, L, L^), being a 
maximum with zero slope, varies only quadratically with deviations from the true L^max, 
and hence may be determined more accurately. Henceforth, the error bars we use for g^ax 
are the above estimated statistical errors, while the error bars we use for max are the above 
defined systematic errors. 

Proceeding in this way at other temperatures, we plot in Fig. |3]the values of gmax(T, L) 
vs. T for L = 8, 10, 12, 16 and 20. The different curves all intersect at a common point, 
determining — 2.05. To determine the correlation length exponent z/, and get a more 
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T 

FIG. 3: Binder ratio maximum gmax(T,L) vs. T for various system sizes L. The common inter- 
section determines Tc — 2.05. Solid lines are linear fits to the data. 

precise estimate of Tc, we now fit the data of Fig. Elto a polynomial expansion as in Eq. fll4p . 
In Table Ewe show the results from both quadratic and cubic polynomial fits, using different 
system sizes L; we systematically drop the smallest sizes since scaling holds only in the 
asymptotic large L limit. Our results give a consistent value of ~ 2.05. The values of u 
that we obtain are consistent within the estimated statistical errors, however we see a small 
systematic increase in the value of u when we restrict the data to larger system sizes. We 
therefore estimate u = 1.2 ±0.1. Note that our value u > 1, satisfies the Chayes lower bound 
condition,— as generalized^ for correlated disorder, u > 2/d*, where d* = 2 is the number of 
dimensions in which the system is disordered. In Fig. 0] we replot the data of Fig. |21 in a 
scaled form, gms^xiT, L) vs. ((T — T^) /Tc)L^^'' . We use the value of and v from TableUfor 
the cubic fit to sizes L = 12 — 20. The solid line is the fitted cubic polynomial. As is seen, 
the data collapse is excellent. 

Having found the value of Tc, we next determine the anisotropy exponent C,. In Fig. Elwe 
show a log-log plot of our data for T^max vs. T, at the temperature T = 2.05 ~ T^. Fitting 
to Eq.jni), 

max 

we get the results summarized in Table |n] for different ranges of 
system sizes L. The results are consistent within the estimated statistical error and we find 
C = 1.3 ± 0.1. To check the consistency of our value for in Fig. IHlwe plot 5'(Tc, L, L^) vs. 
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order 
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8-20 


quadratic 


2.052 ± 0.001 


1.07 ±0.03 




cubic 


2.052 ± 0.001 


1.04 ±0.06 


10-20 


quadratic 


2.052 ± 0.001 


1.10 ±0.04 




cubic 


2.052 ± 0.001 


1.10 ±0.06 


12 - 20 


quadratic 


2.051 ±0.001 


1.18 ±0.04 




cubic 


2.051 ±0.001 


1.18 ±0.05 



TABLE I: Fitting parameters and v from quadratic and cubic polynomial scaling fits to the 
data of Fig. 13 Results for different ranges of system sizes L are shown. 
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-0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 

FIG. 4: Scaling collapse of 5max(r, L) vs. {{J -T^) IT^)L^I'' . The values of = 2.051 and v = 1.18 
from the last row of Tableware used. The solid curve is the fitted cubic polynomial. 

L^/L'', using our data at T = 2.05 and the above determined value of C = 1-3. As expected 
from Eq.(|lip. the data for the different values of L and Lz show a very good collapse to a 
single scahng curve. 

We have also tried a similar scaling analysis for the helicity moduli product, LF'\Tx^ z]-, 
of Eq. ()21|) . However here we have found less satisfactory results. We find that for a given 
system size L, the L^max where L^[Ta;T^] has its maximum occurs at a smaller value of Lz 
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FIG. 5: Log-log plot of L^max vs. L, at T = 2.05 ~ Tc. The solid straight line is the best power 
law fit using the data for L = 10 - 20 and yields the value C = 1-329 ± 0.08 (see TableHH- 
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10 


12 


c 


1.29 ±0.05 


1.33 ±0.08 


1.37±0.12 



TABLE II: Anisotropy exponent ^ from power law fits, Lzmax ~ to system sizes L = Lmin ~ 20. 

than was the case for the Binder ratio Qmax- Such smaller system sizes presumably have 
larger corrections to scaling. We have also found the statistical error in L^[T^T^] to be 
larger than we found for (7max, possibly because the Binder ratio g involves a ratio between 
fluctuating quantities and so has smaller sample to sample fluctuations.-^ As a consequence 
of these two effects, we could not arrive at a convincing determination of Tc and u from 
the L'^lT^Tz] data alone. However, to illustrate our results we can make use of the values 
of Tc ~ 2.05 and ( ~ 1.3 found in our analysis of Qmax- In Fig- 13 we therefore show a 
scaling collapse similar to that of Fig. El plotting L^lT^Tz] vs. L^/L'', using our data at 
T = 2.05 ~ Tc and the above value of (. 

We see clearly in Fig.[2|the above effects: error bars are considerably larger than in Fig. El 
and the peak is at a smaller value of Lz/L^. The scaling collapse is not bad for the bigger 
systems sizes, corresponding to larger values of L^/L^. However it is rather scattered near 
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FIG. 6: Scaling collapse of g{Tc, L, Lz) vs. Lz/L'^, for data at T = 2.05 ~ Tc, using ( = 1.3. The 
solid line is a guide to the eye only. 

the peak and below it. We conclude that it would be necessary to average over many more 
disorder realizations to reduce the errors, and perhaps also go to larger system sizes, in order 
to get a convincing scaling analysis from the helicity product L'^lT^Tz] on its own. 




IV. DYNAMIC CRITICAL BEHAVIOR 



As one approaches the critical temperature T^, we expect relaxation times to diverge as 
T ~ defining the dynamic critical exponent z. To compute equilibrium critical exponents, 
it is sufficient that the simulation dynamics satisfies detailed balance; the details of the 
dynamics are otherwise irrelevant. Thus the exact duality between Tixv and Tiioop implies 
that the unscreened and the strongly screened superconductors have the same equilibrium 
critical behavior. For the dynamic critical behavior, however, the value of z will in general 
depend on the details of the dynamics,^ and some works suggest that it may even vary 
for different types of relaxational dynamics or different boundary conditions . ^''i^^ There is 
thus no reason, a priori, to expect the same dynamic critical behavior for the XY model, 
expressed in terms of a dynamical rule for the phase variables 6i, as compared to the loop 
model, expressed in terms of a dynamical rule for the vortex line variables n^^. In this 
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FIG. 7: Attempted scaling collapse of ^^[T^jTz] vs. LzjLr- ■ Data is for T = 2.05 ~ Tc, using 
C = 1.3. 

section, therefore, we will present results from explicit simulations of the loop model as well 
as the XY model. 

Because the true dynamics of a superconductor is local, it is not physically meaningful 
to compute the dynamic critical exponent within accelerated global algorithms such as the 
Wolff algorithm, which we used to compute equilibrium properties. We therefore will use a 
local Monte Carlo dynamics for both TixY and Tiioop- Even within such local algorithms, it is 
not obvious how universal the dynamical critical behaviors may be. Thus it is unclear that 
our results will correspond to what is seen in experiments. Nevertheless it will be interesting 
to see if the two models give similar or different values of z. 

The relative loss of efficiency that results from using such local algorithms means that 
we will be unable to do as extensive an exploration of the parameter space as we did for our 
equilibrium analysis. But this is not necessary. We can make use of our already obtained 
equilibrium results, and simulate only at the value of T = T^, using system aspect ratios 
Lz = 7L''. For our simulation of TCioop, we will simulate the loop model which is exactly dual 
(see Appendix A, Eq. Mlj) ) to the cosine XY model that we have used in our equilibrium 
simulations, so as to make use of these known values of and (. 
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A. Monte Carlo Methods and Scaling 

For the XY model of an unscreened superconductor we use a standard single spin heat 
bath algorithm, with fixed periodic boundary conditions on the di. In this algorithm, a phase 
angle 9i is selected at random and replaced with a new randomly chosen 9[. This update 
attempt is then accepted with probability + exp(A£'/T)] where /S.E is the change in 
energy. One sweep, consisting of = L'^L^ update attempts, is taken as one time step. 
At = 1. We average over 300 — 700 disorder realizations depending on system size L. 

For the loop model of a strongly screened superconductor, we again use a heat bath 
algorithm in which the attempted excitation consists of an elementary vortex loop circulating 
about a randomly chosen plaquette of the grid. Adding only such closed loop excitations 
corresponds to the ensemble in which the average internal magnetic field is constrained to 
B = (see Appendix A). One sweep, consisting of 3A^ such update attempts, is taken as one 
time step. At = 1. We average over 1000 — 2000 disorder realizations depending on system 
size L. 

In general, we expect the relaxation time r to obey the scaling equation, 

r(T, L, L,) = (t6l/^ Ljh^) , (22) 

where h is an arbitrary length rescaling factor. For b = L, T = Tc, and Lz = 'jL'', this 
reduces to the simple, 

Tr^ . (23) 

For both the XY model and the loop model, we simulate with values of Lz = as 
determined by the fit shown in Fig. For the XY model, to approximate non-integer values 
of Lz we use linear interpolation of simulation data for the two closest integer values of Lz- 
For the loop model we simply use results from the closest integer value of Lz- 

B. Observables 

L XY Model 

For the XY model we have tried two independent methods of determining z, analogous to 
the two quantities g and L^[Ta.T^] used in our equilibrium simulations. The first is to look 
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at the decay of correlations in the order parameter M of Eq. ()16|) . defining the relaxation 
time r by, 



to 



t=i 



{M{t)M{0)) 



(24) 



where to is chosen large enough so that r is independent of to- The ratio in the above ensures 
that the quantity being summed over has scaling dimension zero, and hence the sum scales 
as r ~ L^. 

The second method is to look at correlations of the supercurrent defined by, 
. dUxY 



OA, 



Aa = 



sm 



(25) 



In terms of one can define the conductance in the fi direction by the Kubo formula,-^ 



(26) 



t=-to 



where again to is chosen large enough that is independent of to- Since = 
{dTCxY /dAfj)\^^^Q, and Hxy and are scale invariant, then J^, and hence the corre- 
lation summed over in in the definition of G^, has scaling dimension zero. Therefore, the 



sum which defines G/j, scales as r ~ L"^. 



2. Loop Model 

For the loop model we consider the total resistance, defined as follows?^ Let Qfi{t) be 
the total projected loop area with normal in direction /t at simulation time t. Each time an 
oriented elementary vortex loop with normal in direction ±/i is accepted, changes by ±1. 
Let AQ^it) = Q^(t) — Q^{t — 1) be the total change in this area after one sweep through 
the entire system; each sweep represents At = 1. In one such sweep, the total average phase 
angle change across the length of the system (in the dual screened XY superconductor model) 
in direction (1 is just 2'KAQ^/LyL„^ where /x, z/, a are a cyclic permutation of x, y, z. By 
the Josephson relation, the total voltage drop across the system in direction fi will then be 
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^min 


6 


8 


10 


12 


ZXY 


2.72 ±0.04 


2.69 ±0.05 


2.63 ± 0.07 


2.60 ±0.03 



TABLE III: Dynamic exponent zxy from power law fits, t ^ , to system sizes L = Lmin ~ 20. 

Henceforth we define our units of voltage such that h/2e = 1. We then define the total 
resistance in direction fi by the Kubo ioTumla^ 

= ^ E mv,m,m] , (28) 

t=-to 

where again to is chosen large enough so that i?^ is independent of to. Since the total 
voltage drop is the time rate of change of the total phase angle difference across the 
system, and since the total phase angle difference is a scale invariant quantity, we have the 
scaling ~ 1/r. Thus the resistance above scales as, 

~ 1/r ~ L-^ . (29) 

C. Results 

1. XY Model 

In Fig. IHl we show a log-log plot of our results for the order parameter relaxation time r 
of Eq. ()24|) vs. system size L, for T = 2.05 ~ and ~ L^. Our results are obtained using 
5 X 10^ MC sweeps to equilibrate, followed by 10® sweeps to compute averages. Fitting to 
the power law, r ~ L^, we get the results summarized in Table ITTTl for different ranges of 
systems size L. The results are consistent within the estimated statistical error, with a small 
systematic tendency to lower values as we restrict the fitted data to larger system sizes. We 
find Zxy = 2.63 ± 0.07. 

As another check on our above determination of zxy, we consider the following. In 
principal, r is defined by taking tg in Eq.(j2D) sufficiently large so that r is independent of 
to; our data in Fig. |H1 satisfies this condition. How big tg must be for this to happen is set 
by the time scale r itself. Therefore, we expect that if we compute r for arbitrary to, then 
r(to) should scale as, 

r(to) ~ L^f (to/r) ~ L^f (to/L^) . (30) 
19 




FIG. 8: Log-log plot of order parameter relaxation time r of Eg . (124(1 vs. system size L, for 
T = 2.05 ~ Tc and ~ L''. Solid line is the best power law fit for sizes L = 10 — 20, and 
determines z = 2.63 ± 0.07 (see Table Iffll). 

In Fig. El we show a log-log plot of T{to)/L^ vs. to/L^ for various sizes L (again using 
T = 2.05 ~ Tc and ~ L'^). Choosing the value zxv = 2.63 obtained from the fit in Fig. |H1 
we find an excellent collapse of all the data. For large t^/L^ we see that the curve does 
indeed saturate to a finite constant as expected, however the collapse holds for the entire 
range of to- 

Finally, we plot in Fig.lTUlthe conductances of Eq.ipm). andG^ vs. L, forT = 2.05 ^ Tc 
and Lz ~ L^. Our results are for 2 x 10^ MC sweeps to equilibrate, followed by 4 x 10^ sweeps 
to compute averages. Fitting to the power law, ~ L^f*, we get the results summarized 
in Table IIVI for different ranges of systems size L. For G^ the results z ~ 2.66 ± 0.04 
are consistent, within errors, with that obtained from our analysis of the order parameter 
relaxation time r. For Gz, we get values for z that are somewhat larger. However if one 
compares the data points for Gx and Gz directly, one sees that the values are all roughly 
equal within the estimated error, except for the smallest size L = 4 (probably too small to 
be in the scaling limit) and for the largest size L = 20. Our fit for z from the Gz data is 
skewed by this one L = 20 data point. If we restrict our fit to sizes L = 8 — 16, we then find 
Zz = 2.82 ± 0.03. This is still somewhat larger than what we get from G^, but within two 
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FIG. 9: Log-log scaling plot of order parameter relaxation time T{tQ)/L^ vs. tg/L^ for T = 2.05 ~ 
Tc, Lz ~ -L^, and various values of L. Using zxy = 2.63 gives an excellent collapse for the entire 
range of to- 
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8 


10 


12 




2.71 ±0.02 


2.75 ±0.03 


2.66 ± 0.04 


2.69 ±0.06 


Zz 


2.77 ±0.02 


2.87 ±0.03 


2.91 ± 0.04 


3.06 ± 0.07 



TABLE IV: Dynamic exponent zxy from power law fits, ~ L^^^, to system sizes L = Lmin — 20. 
standard deviations of for the same range of sizes L = 8 — 16. 



2. Loop Model 

For our foop simulations we use the interaction of Eq. pT|) . exactly dual to our XY model. 
This interaction is computed using the same distribution of Jj^ as we used for the XY model, 
and we simulate at the same value of T = 2.05 as gives the critical point of the XY model. 
We also use the same values of = 7-^^ as we used for the XY model, as determined from 
Fig. El In Fig. ^2 we give our results for the resistance of the loop model, Eq. ()28|) . as a 
log-log plot of Rx and Rz vs. system size L. Our results are from 12 x 10^ MC sweeps 
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FIG. 10: Log-log plot of conductances Gx and vs. L, for T = 2.05 ~ Tc and ~ L'-. The 
fitted straight hnes determine Zx = 2.66 and = 2.91. 

to equilibrate, followed by 24 x 10^ sweeps to compute averages. Fitting to the power law 
of Eq.()29|). Rx ~ "we get the results summarized in Table IV| for different ranges of 
systems size L. The results are consistent within the estimated statistical error, and we find 
zioop = 3.4±0.1. 

For the case of Rz-, parallel to the columnar defects, our simulations were not sufficiently 
long to observe the necessary saturation of Rz{to) with increasing to, except for the smallest 
system sizes L < 12. We do not believe that any estimate of 2;ioop based on such small 
system sizes would be meaningful. We can, however, perform the following consistency 
check. Similar to our discussion concerning r(to) (see Eq. ljHUj) ). we can compute R^ of 
Eq . (j28|) for finite times t^, and we expect R^{to)L^ to scale with the variable to/L^. In 
Fig. ^1 we make such a log-log scaling plot using the value of z = 3.4 found for Rx in 
Fig. ^2 For Rx we see that the collapse is excellent for all times to 5 and the scaling curve 
saturates to a constant at large to/L^ as expected. For Rz, we find a good collapse for all 
but the largest times. We see that Rz{to) saturates only for the smallest systems, and it 
is only here that the collapse appears to be breaking down. We conclude that these small 
system sizes are not large enough to expect scaling for Rz to hold. 

We can also try to independently determine the dynamic exponent z by fitting to a data 
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^min 


6 


8 


10 


12 


^loop 


3.23 ±0.05 


3.33 ± 0.07 


3.39 ±0.11 


3.38 ±0.14 



TABLE V: Dynamic exponent zioop from power law fits, ^ L ^, to system sizes L = Lmin ~ 20. 

collapse as in Fig. ^Jfor all times to? rather than just the asymptotic large time limit. The 
inset to Fig. El shows the resulting of such fits as the fitting parameter z is varied. For 
Rxi the shows a sharp minimum at z = 3.45, in good agreement with our earlier value 
of 2; = 3.4 from Fig. ^2 For R^, the has a minimum at the somewhat higher value of 
z = 3.7, however the minimum is very shallow, indicating a relative insensitivity of the data 
to variations in z. We conclude that both the data for i?^ and Rz are consistent with a 
dynamic exponent 2;ioop = 3.4 ± 0.1. 

10-3 
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N 10-5 
pit 10-6 
10-7 
10-8 

4 5 6 7 8 9 10 ^ 20 30 

FIG. 11: Log-log plot of resistance Rx and Rz of the loop model vs. L. The solid line is the best 
power law fit, Rx ~ L~^, for sizes L = 10 — 20, and determines the value zioop = 3.4 ±0.1 (see 
Table EI). 
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FIG. 12: Log-log scaling plot of time dependent resistance Rx{to)L'^ and Rz{tQ)L^ of the loop 
model vs. to/L^. The value z = 3.4 obtained from the fit in Fig. ^Jis used. The inset gives the 
error of the data collapse, as the exponent z is varied. 

V. DISCUSSION AND CONCLUSIONS 



We have studied the equilibrium and dynamic critical behavior of the zero magnetic field 
superconducting phase transition for a type-II superconductor with quenched columnar dis- 
order. We have considered both the "unscreened" XY model in which Aq — > oo, and the 
"strongly screened" loop model in which Aq ~ ^q. A duality transformation establishes that 
these two models are in the same equilibrium universality class. Using numerical simula- 
tions of the XY model, we find, in agreement with a generalized Harris criterion, that the 
universality class of the transition is different from the pure model, and we find that scaling 
is anisotropic. We find the value for the correlation length exponent, u = 1.2 ± 0.1, and for 
the anisotropy exponent, C = 1.3 ± 0.1. 

Using the value of the critical temperature and the anisotropic scaling determined from 
the equilibrium analysis, we carry out simulations at the critical point to determine the 
dynamic critical exponent z of both the XY and loop models for local Monte Carlo dynamic 
rules. For the "unscreened" XY model, with a single spin heat bath dynamics, we find 
zxY = 2.6 ±0.1. For the "strongly screened" loop model, with a heat bath dynamics applied 
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to elementary loop excitations, we find 2;ioop = 3.4 ± 0.1. 

A similar random 3D XY model has been studied by Cha and Girvin^ in the context of 
the quantum phase transition in the two dimensional boson Hubbard model. In their model 
disorder was introduced as uniformly distributed random bonds in the z (imaginary time) 
direction, Jj^, so as to model bosons with random charging energy. They found equilibrium 
critical exponents z/ = 1.0 ± 0.3 and C = 1-07 ± 0.03 (our anisotropy exponent C for the 
classical 3D model is equivalent to their "quantum dynamic exponent" z for the 2D quantum 
problem). However their analysis for such a system with anisotropic scaling, C > !> was based 
on a more ad hoc approach of (i) trying various values of C and seeing which appeared to give 
the best data collapse for systems of different size L, and (ii) measuring real space correlations 
in a system of a fixed size and fitting to assumed power law decays. Their largest system 
size, 16^ X 15 is also smaller than ours and they do not use the Wolff algorithm to accelerate 
their equilibration. While it is possible that introducing the randomness differently (along 
z rather than in the xy plane) might effect the universality class, we believe it is more likely 
that this is not the case, and that our results are more systematic and hence more accurate 
than those of Cha and Girvin. 

Prokof'ev and Svistunov^° have simulated the loop model of Eq.(jS)) in the context of the 
same two dimensional disordered boson problem as Cha and Girvin. For their "off- diagonal" 
disorder case they put the disorder into the bonds along the z direction, making their model 
dual to that of Cha and Girvin. They report an anisotropy exponent C = 1.5 ± 0.2, which 
agrees with ours within the estimated errors. They were unable to determine the correlation 
length exponent v. We note that while they use an accelerated "worm" algorithm and have 
good statistics for quite large system sizes, they determine their exponents by fitting to 
real space correlation functions for their biggest size system, as did Cha and Girvin, rather 
than doing any systematic finite size scaling that takes into account the anisotropic scaling 
present in the model. 

Experimental investigation of the zero field transition with colummnar disorder has been 
undertaken by Kotzler and co-workera^ii for YBaCuO thin films. Measuring the frequency 
dependent conductivity transverse to the columnar disorder, which is expected to scale asi^ 
ax(cu,r) = t^'^-^>d{ujt'"') (where t = {T - T^)/T^), they find^ the combinations uC = 1.7 
and z/( = 5.53. This compares with our values = 1.56, and z/( = 2 for the unscreened 
XY model, and z/( = 2.6 for the strongly screened loop model. Our value of is conceivably 
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consistent with the experimental value, within possible errors. However both of our values 
of zjC, seem too small. It may be that our simple local Monte Carlo dynamics does not 
adequately capture the true dynamics of a real superconductor. On the other hand, if we 
use our value of C = 1-3, then Kotzler's results imply a dynamic exponent oi z = 7.2, which 
seems extraordinarily large. 

We may also compare our dynamic exponents with those obtained from the disorder free 
model. For the the strongly screened limit of the loop model, Lidmar et al^ find the value 
z\oo^ ~ 2.7; moreover they find this value to be insensitive to the presence of uncorrelated 
point disorder. For relaxational dynamics of the phase angle variable in the XY model, a 
value of 2; ~ 2 is expected,— and this is what was found in numerical simulations by Jensen 
et al.'^'^ using a method similar to our scaling of conductance, Eq.(j2ni)- The result 2;ioop > zxy 
thus seems common for both the pure and columnar disordered cases. 

In our work we have considered only simple relaxational dynamics for the phase angles of 
the unscreened XY model. Two other possible dynamics might be considered. One would 
be to do a loop dynamics, similar to what we have done here for the strongly screened loop 
model, only now as applied to the strongly interacting loops of the unscreened XY model. 
The other would be to use resistively shunted junction (RSJ) dynamics for the phase angles 
of the XY model. Both such approaches have been previously used for the disorder free 
case. For both loop dynamicsi^*^ and RSJ dynamic o^^i^^ the dynamic exponent z ~ 1.5 
was found, smaller than the value obtained by simple phase angle relaxational dynamics. 
Investigating these other dynamics for the case of columnar disorder remains for future work. 
We only note here that if the above trend remains true for columnar disorder, and that these 
other dynamics reduce z from that of relaxational dynamics, then it becomes even harder 
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to explain the large value of z( observed experimentally in Ref. 

Finally, we note that similar equilibrium exponents to those found in this work were 
also found for the case of an unscreened superconductor with columnar defects in a finite 
applied magnetic field. For that case the values^ u = 1.0 ± 0.1 and ( = 1.25 ± 0.1 were 
found. Although these are close to the values we find here for zero applied field, there is 
no apparent reason that the zero and finite field cases should be in the same universality 
class. We also note that once a finite field is applied, the duality between the unscreened 
and strongly screened superconductor models, that exists for zero field, breaks down. 
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Appendix A 

In this section we review the duality transformatio n^i^^i^^ from ?ixY of Eq. (Q) to TYioop of 
Eq.©. Consider first a general 27r periodic interaction Vi^{(f)) instead of the —{Ji^/T) cos(0) 
of Eq.(0). For the generalized fixed twist boundary condition and the corresponding Hamil- 
tonian of Eq.(j3]), we can write the partition function as, 

Z = (lll'^ e-^-^-(^^''^+M-^MAM) . (31) 
where the 6'^ obey periodic boundary conditions. Defining the Fourier transform Vj^ by, 

oo 



n,-„=— oo 



and substituting into Eq. ljHTj) gives 

"Stt 7/1/ 



= ^ e'^^M[^.MKM)+-..A,/L,] I ^'"^ j ^ _ (34) 

One is now free to do the integrals over the 6j. The result is a product of Kronecker deltas 
constraining the variables rij^ to be divergenceless, as in Eq.©. Defining the "winding 
numbers" by, 

^M^T^E^^A^' (35) 



we get, 

Z= ^'e-^^M^^MKM)-iE^H/^A^ ^ (3g) 
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where the prime on the summation denotes the divergenceless constraint of Eq.( 
A common choice for Vifj,{(f)) is the Villain interaction,— 



oo 



J, 



e-v.^W = J2 e-^('^-2"'")' . (37) 

m=— oo 

In this case one has for its transform, 

Ki^) = T^ri' ■ (38) 

The partition function of Eq. ()36|) . with periodic boundary conditions = 0, then becomes, 

Z= ^ V^^^'^^^'^^> , (39) 

with 

g,,/f = T/J,, . (40) 

The above is just a model of short ranged interacting loops with onsite repulsion ~ and 
inverted temperature scale T ~ 1/T. 

For our simulatons, with Vi^(0) = —{Jifj,/T) cos(0), one hast^^ 

^ /„(j.^/T) , (41) 

where /n(x) is the modified Bessel function of the first kind. Since /„(x) is an increasing 
function of \n\ for fixed x, the above similarly gives a short ranged loop model with onsite 
repulsion. It is this interaction of Eq.()4H) that we use in our dynamic simulations of the 
loop model in Section IV. 

We can now demonstrate several interesting results concerning phase coherence in the 
XY model, by considering the behavior as a function of the twist A^. The XY model is 
phase coherent when the total free energy JF varies with A^. Using JF(A^) = — TlnZ(A^), 
and Eq. ljHUj) above, we find, 

^{W,)o , (42) 



where (. . .)o indicates an average in the ensemble with A^ = 0. Now since dJ-'/dAfj_ must be 
a real quantity (as may be seen by considering its evaluation in the original XY model Tixv 
of Eq.(jH)), and since (W^)o must similarly be real (as may be seen by considering Tiioop), 
the only way for Eq. ljl^ to hold is if 9JF/9A^|^ = (VF^)o = 0. This then demonstrates 
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that = 0, i.e. periodic boundary conditions on the 9i, is the twist that minimizes the 
free energy. 

Finally, returning to Eq.(|36p. we note that in the fluctuating twist ensemble^ for the 
XY model, in which is averaged over as a thermally fluctuating degree of freedom, 
the corresponding loop model obeys the additional constraint of zero winding, VF^ = 0, in 
each individual configuration. When viewing Tiioop as the Hamiltonian of vortex loops in 
a strongly screened superconductor, this corresponds to the ensemble in which the average 
internal magnetic field is constrained to vanish, 5^ = 0, in each configuration. 
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